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ABSTRACT 

This  paper  is  focused  on  the  dynamic  formulation  of  mechanical  joints  using  different 
approaches  that  lead  to  different  models  with  different  numbers  of  degrees  of  freedom.  Some  of 
these  formulations  cdlow  for  capturing  the  joint  deformations  using  discrete  elastic  model  while 
the  others  are  continuum-based  and  capture  joint  deformation  modes  that  cannot  be  captured 
using  the  discrete  elastic  joint  models.  Specificcdly,  three  types  of  joint  formulations  are 
considered  in  this  investigation;  the  ideal,  compliant  discrete  element,  and  compliant  continuum- 
based  joint  models.  The  idecd  joint  formulation,  which  does  not  cdlow  for  deformation  degrees  of 
freedom  in  the  case  of  rigid  body  or  small  deformation  analysis,  requires  introducing  a  set  of 
algebraic  constraint  equations  that  can  be  handled  in  computational  multibody  system  (MBS) 
algorithms  using  two  fundamentally  different  approaches:  constrained  dynamics  approach  and 
penalty  method.  When  the  constrained  dynamics  approach  is  used  the  constraint  equations  must 
be  satisfied  at  the  position,  velocity,  and  acceleration  levels.  The  penalty  method,  on  the  other 
hand,  cannot  ensure  that  the  cdgebrcdc  equations  are  satisfied  at  the  acceleration  level.  In  the 
compliant  discrete  element  joint  formulation,  no  constraint  conditions  are  used;  instead  the 
connectivity  conditions  between  bodies  are  enforced  using  forces  that  can  be  defined  in  their 
most  general  form  in  MBS  algorithms  using  bushing  elements  that  allow  for  the  definition  of 
genercd  nonlinear  forces  and  moments.  The  new  compliant  continuum-based  joint  formulation, 
which  is  based  on  the  finite  element  (FE)  absolute  nodal  coordinate  formulation  (ANCF),  has 
severed  advantages:  (I)  It  captures  modes  of  joint  deformations  that  cannot  be  captured  using 
the  compliant  discrete  joint  models;  (2)  It  leads  to  linear  connectivity  conditions,  thereby 
allowing  for  the  elimination  of  the  dependent  variables  at  a  preprocessing  stage;  (3)  It  leads  to  a 
constant  inertia  matrix  in  the  case  of  chain  like  structure;  and  (4)  It  automatically  captures  the 
deformation  of  the  bodies  using  distributed  inertia  and  elasticity.  The  formulations  of  these  three 
different  joint  models  are  compared  in  order  to  shed  light  on  the  fundamental  differences 
between  them.  Numerical  results  of  a  detailed  tracked  vehicle  model  are  presented  in  order  to 
demonstrate  the  implementation  of  some  of  the  formulations  discussed  in  this  investigation. 
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PROBLEM  DEFINITION 


Ideal  Joint  Formulation 


Accurate  formulation  of  mechanical  joints  is 
necessary  in  the  computer  simulation  of 
multibody  system  (MBS)  that  represent 
many  technological  and  industrial 
applications.  An  example  of  these  MBS 
applications,  in  which  accurate  modeling  of 
joint  compliance  is  necessary,  is  the  tracked 
vehicle  shown  in  figure  1.  The  links  of  the 
track  chains  of  this  vehicle  are  connected  by 
pin  joints  that  can  be  subjected  to  significant 
stresses  during  the  vehicle  functional 
operations.  Nonetheless,  there  are  different 
joint  formulations  that  can  lead  to  different 
dynamic  models  which  have  different 
numbers  of  degrees  of  freedom.  This  study 
investigates  the  use  of  three  different 
methods  for  formulating  mechanical  joints 
in  MBS  applications.  These  three  methods 
are  the  ideal,  the  compliant  discrete  element, 
and  the  compliant  continuum-based  joint 
formulations.  These  three  different  methods 
are  described  below. 


Chassis 


Sprccket 


Road-wheel 


Figure  1:  Ml  13  tracked  vehicle  model  in 
SAMS/2000 


The  ideal  joint  formulation  is  based  on  a  set 
of  algebraic  equations  that  do  not  account 
for  the  joint  flexibility;  this  is  regardless  of 
whether  or  not  the  body  is  flexible.  The 
algebraic  joint  equations  are  expressed  in 
terms  of  the  coordinates  of  the  two  bodies 
connected  by  the  joint.  These  algebraic 
equations  are  considered  as  constraint 
equations  which  can  be  enforced  using  two 
fundamentally  different  methods;  the 
constrained  dynamics  approach  and  the 
penalty  method.  In  the  constrained  dynamics 
approach,  the  technique  of  Lagrange 
multipliers  or  a  recursive  method  is  used.  In 
this  case,  the  joint  constraint  equations  must 
be  satisfied  at  the  position,  velocity,  and 
acceleration  levels.  The  number  of  degrees 
of  freedom  of  the  model  in  this  case  is  equal 
to  the  number  of  the  system  coordinates 
minus  the  number  of  the  algebraic  joint 
constraint  equations.  In  the  penalty  method, 
on  the  other  hand,  the  number  of  degrees  of 
freedom  of  the  model  is  not  affected  by  the 
number  of  joint  constraint  equations.  These 
joint  algebraic  equations  are  enforced  using 
high  stiffness  penalty  coefficients  that 
ensure  that  the  algebraic  constraint 
equations  are  satisfied  at  the  position  level. 
The  penalty  method  does  not  ensure  that  the 
constraint  equations  are  satisfied  at  the 
acceleration  level. 
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Compliant  Discrete  Element  Joint 
Formulation 

In  this  approach,  no  algebraic  equations  are 
used  to  describe  the  joints  between  bodies  in 
the  system.  The  connectivity  between  bodies 
is  described  using  force  elements  that  have 
forms  defined  by  the  user  of  the  MBS  code. 
MBS  system  codes  have  bushing  elements 
that  can  be  used  to  define  general  linear  or 
nonlinear  force  and  moment  expressions. 
The  stiffness  and  damping  coefficients  in  the 
force  and  moment  expressions  can  be 
selected  by  the  user.  The  busing  elements 
can  be  used  to  model  the  joint  compliance  in 
the  case  of  rigid  and  flexible  body  dynamics. 
It  is  important,  however,  to  point  out  that 
adding  bushing  elements  has  no  effect  on 
the  number  of  degrees  of  freedom  of  the 
model.  Unlike  the  penalty  method,  the  use 
of  bushing  element  does  not  require  the 
formulation  of  algebraic  joint  equations. 
Bushing  elements  allow  for  systematically 
introducing  three  force  components  and 
three  moment  components. 


Compliant  Continuum-Based  Joint 
Formulation 

The  new  finite  element  (FE)  absolute  nodal 
coordinate  formulation  (ANCF)  allows  for 
systematically  developing  new  joint 
formulations  that  capture  modes  of 
deformation  that  cannot  be  captured  using 
the  discrete  joint  models.  It  also  allows  for 
modeling  body  flexibility  using  new  FE 
meshes  that  have  constant  inertia  and  linear 
connectivity  conditions.  Specifically,  the 
compliant  ANCF  continuum-based  joint 
formulation  has  the  following  advantages: 

1.  ANCF  finite  elements  allow  for 
developing  new  joint  formulations  that 


capture  deformation  modes  that  cannot  be 
captured  using  compliant  discrete  joint 
formulations.  The  use  of  the  ANCF 
gradient  coordinates  allows  for 
developing  different  joint  models  with 
different  numbers  of  degrees  of  freedom 
that  allow  for  different  strain  modes. 

2.  The  use  of  the  ANCF  gradient  coordinates 
allows  for  developing  linear  joint 
constraint  equations.  These  linear 
algebraic  equations  can  be  used  to 
eliminate  dependent  variables  at  a 
preprocessing  stage,  thereby  significantly 
reducing  the  model  dimensionality. 

3.  ANCF  finite  elements  can  also  be  used  to 
model  the  body  deformation  in  addition 
to  the  joint  compliance.  Distributed 
inertia  and  elasticity  are  used  for  both 
body  flexibility  and  joint  compliance. 

4.  ANCF  finite  elements  lead  to  new  types 
of  FE  meshes  that  have  constant  inertia,  a 
feature  that  can  be  exploited  to  develop  a 
sparse  matrix  structure  of  the  MBS 
dynamic  equations. 

It  is  the  objective  of  this  investigation  to 
provide  a  comprehensive  study  of  different 
joint  formulations  and  demonstrate  the 
fundamental  differences  between  them  when 
applied  to  the  analysis  of  complex  tracked 
vehicle  system  models.  Simulation  models 
of  the  Ml  13  tracked  vehicle  will  be  used  to 
compare  the  numerical  results  obtained 
using  the  joint  formulations  discussed  in  this 
paper.  These  results  are  obtained  using  the 
general  purpose  MBS  computer  code 
SAMS/2000  [17]  which  allows  for 

systematically  modeling  MBS  applications 
using  the  augmented  formulation,  penalty 
method,  bushing  elements,  and  ANCF  finite 
elements.  In  the  augmented  formulation,  the 
technique  of  Fagrange  multipliers  is  used. 
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The  computational  algorithm  used  in 
SAMS/2000  ensures  that  the  algebraic 
constraint  equations  are  satisfied  at  the 
position,  velocity,  and  acceleration  levels. 


TRACKED  VEHICLES: 

BACKGROUND 

High  mobility  tracked  vehicle  such  as 
military  battle  tanks  and  armored  personal 
carriers  are  designed  for  the  mobility  over 
rough  and  off-road  terrains.  Investigations 
on  the  dynamic  analysis  of  such  tracked 
vehicles  shown  in  figure  1  have  been  limited 
because  of  the  complexity  of  the  forces 
resulting  from  interaction  between  the 
vehicle  components.  These  forces  are 
impulsive  in  nature,  and  their  dynamic 
modeling  requires  sophisticated 

computational  capabilities.  Several  two 
dimensional  models  for  the  analysis  of 
tracked  vehicle  have  been  developed. 
Galitsis  [3]  demonstrated  that  the  dynamic 
track  tension  and  suspension  loads  in  high 
speed  tracked  vehicles  developed  by 
analytical  methods  are  useful  in  evaluating 
the  dynamic  characteristics  of  the  tracked 
vehicle  components.  Bando  et  al.  [1] 
outlined  a  procedure  for  the  design  and 
analysis  of  rubber  tracked  small-size 
bulldozers,  and  presented  a  computer 
simulation  which  was  used  in  the  evaluation 
of  the  vehicle  performance.  Both  steel  and 
continuous  rubber  tracks  are  modeled  by 
discretizing  them  into  several  rigid  bodies 
connected  by  compliant  elements.  The 
simulation  results  indicate  that  the  vehicle 
has  favorable  characteristics,  such  as  less 
damage  to  the  road  surface,  and  reduced 
vibration  and  noise.  Murray  and  Canfield  [7] 
used  general  purpose  multibody  computer 
codes  to  model  a  simple  track  link  and 
sprocket  system.  The  behavior  of  the 


interaction  between  the  track  link  and  the 
sprocket  was  illustrated  graphically  and  it 
was  found  that  the  computer  time  can  be 
significantly  reduced  by  using 
supercomputers.  Nakanishi  et  al.  [8] 
developed  a  two  dimensional  contact  force 
model  for  planar  analysis  of  multibody 
tracked  vehicle  systems.  Modal  parameters 
(modal  mass,  modal  stiffness,  modal 
damping,  and  mode  shapes),  which  are 
determined  experimentally,  are  employed  to 
simulate  the  nonlinear  dynamic  behavior  of 
a  multibody  tracked  vehicle  which  consists 
of  interconnected  rigid  and  flexible 
components.  The  equations  of  motion  of  the 
vehicle  are  formulated  in  terms  of  a  set  of 
modal  and  reference  generalized 
coordinates,  and  the  theoretical  basis  for 
extracting  the  component  modal  parameters 
of  the  chassis  from  the  modal  parameters  of 
the  assembled  vehicle  is  described. 

A  number  of  approaches  have  been 
proposed  in  the  literature  for  developing 
three-dimensional  MBS  models.  Choi  et  al. 
[2]  developed  the  nonlinear  dynamic 
equations  of  motion  of  the  three- 
dimensional  multibody  tracked  vehicle 
systems,  taking  into  consideration  the 
degrees  of  freedom  of  the  track  chains.  To 
avoid  the  solution  of  a  system  of  differential 
and  algebraic  equations,  the  recursive 
kinematic  equations  of  the  vehicle  are 
expressed  in  terms  of  the  independent  joint 
coordinates.  In  order  to  take  advantage  of 
sparse  matrix  algorithms,  the  independent 
differential  equations  of  the  three- 
dimensional  tracked  vehicles  are  obtained 
using  the  velocity  transformation  method. 
Three-dimensional  nonlinear  contact  force 
models  that  describe  the  interaction  between 
the  track  links  and  the  vehicle  components 
such  as  the  rollers,  sprockets,  and  idlers  as 
well  as  the  interaction  between  the  track 
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links  and  the  ground  are  developed  and  used 
to  define  the  generalized  contact  forces 
associated  with  the  vehicle  generalized 
coordinates.  A  computer  simulation  of  a 
tracked  vehicle  in  which  the  track  is 
assumed  to  consist  of  track  links  connected 
by  a  single  degree  of  freedom  revolute  joint 
is  presented  in  order  to  demonstrate  the  use 
of  the  formulations  presented  in  their  study. 
Ryu  et  al.  [13]  developed  compliant  track 
link  models  and  investigated  the  use  of  these 
models  in  the  dynamic  analysis  of  high¬ 
speed,  high-mobility  tracked  vehicles.  The 
characteristics  of  the  compliant  elements 
used  in  this  investigation  to  describe  the 
track  joints  are  measured  experimentally.  A 
numerical  integration  method  having  a 
relatively  large  stability  region  is  employed 
in  order  to  maintain  the  solution  accuracy, 
and  a  variable  step  size  integration  algorithm 
is  used  in  order  to  improve  the  efficiency. 
The  dimensionality  problem  is  solved  by 
decoupling  the  equations  of  motion  of  the 
chassis  and  track  subsystems.  Recursive 
methods  are  used  to  obtain  a  minimum  set  of 
equations  for  the  chassis  subsystem.  Several 
simulations  scenarios  including  an 
accelerated  motion,  high-speed  motion, 
braking,  and  turning  motion  of  the  high- 
mobility  vehicle  are  tested  in  order  to 
demonstrate  the  effectiveness  and  validity  of 
the  methods  proposed.  Ozaki  and  Shabana 
[9-10]  evaluated  the  performance  of 
different  formulations  using  a  tracked 
vehicle  model  that  is  subjected  to  impulsive 
forces.  They  developed  models  for  joint 
constraints  and  the  impulsive  contact  forces 
that  result  from  the  interaction  between  the 
track  chains  and  the  vehicle  components  as 
well  as  the  interaction  between  these  chains 
and  the  ground.  The  nonlinear  contact  force 
models  used  in  their  numerical  study  were 
developed,  and  the  formulations  of  the 
generalized  forces  associated  with  the 


generalized  coordinates  used  in  each  of  the 
formulations  were  presented.  Ryu  et  al  [14] 
investigated  the  nonlinear  dynamic 
modeling  methods  for  the  virtual  design  of 
tracked  vehicles  by  using  MBS  dynamic 
simulation  techniques.  The  results  include 
high  oscillatory  signals  resulting  from  the 
impulsive  contact  forces  and  the  use  of  stiff 
compliant  elements  to  represent  the  joints 
between  the  track  links.  Each  track  link  is 
modeled  as  a  body  which  has  six  degrees  of 
freedom,  and  two  compliant  bushing 
elements  are  used  to  connect  track  links. 
Efficient  contact  search  and  kinematics 
algorithms  in  the  context  of  the  compliance 
contact  model  are  developed  to  detect  the 
interactions  between  track  links,  rollers, 
sprockets,  and  ground  for  the  sake  of  speedy 
and  robust  solutions.  Rubinstein  and  Hitron 
[12]  developed  a  three-dimensional  multi¬ 
body  simulation  model  for  simulating  the 
dynamic  behavior  of  tracked  off-road 
vehicles  using  the  LMS-DADS  simulation 
program.  Each  track  link  is  considered  a 
rigid  body  and  is  connected  to  its 
neighboring  track  link  via  a  revolute  joint. 
The  road-wheel  track-link  interaction  is 
described  using  three-dimensional  contact 
force  elements,  and  the  track-link  terrain 
interaction  is  modeled  using  a  pressure - 
sinkage  relationship. 


SCOPE  AND  OBJECTIVES  OF  THIS 
INVESTIGATION 

While  some  of  the  formulations  used  in  the 
analysis  of  tracked  vehicles  require  the  use 
of  MBS  algorithms  that  are  designed  for 
solving  systems  of  differential/algebraic 
equations  arising  from  kinematic  joint 
constraints,  other  formulations  do  not 
require  the  use  of  such  algorithms  but  use 
penalty  forces  or  complaint  elements  to 
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model  the  chain  dynamics.  The  track  links 
are  subjected  to  high  contact  forces  that 
produce  high  stress  levels,  which  often 
cause  damage  to  the  track  links  during  the 
functional  operation  of  the  vehicle.  The 
objective  of  this  paper  is  to  present  and 
compare  between  different  joint 
formulations  that  can  be  used  in  the  dynamic 
modeling  of  complex  tracked  vehicle 
systems.  Better  understanding  of  these 
formulations  can  lead  to  more  accurate,  and 
possibly  faster,  computer  simulations  that 
can  be  the  basis  for  more  reliable 
performance  evaluation  of  the  vehicles.  The 
tracked  vehicles  considered  in  this 
investigation  are  assumed  to  consist  of 
interconnected  bodies  that  can  have  arbitrary 
displacements.  In  the  first  chain  formulation, 
referred  to  in  this  paper  as  the  ideal  joint 
chain  model ,  kinematic  joints  between  the 
track  links  are  described  using  nonlinear 
constraint  equations  that  lead  to  significant 
reduction  in  the  number  of  vehicle  degrees 
of  freedom.  This  joint  model  does  not 
require  assuming  stiffness  and  damping  for 
the  track  link  connectivity,  and  therefore, 
does  not  allow  for  flexibility  between  the 
tack  links;  it  requires,  however,  the  solution 
of  a  system  of  differential  and  algebraic 
equations  if  redundant  coordinate 

formulations  are  used.  Redundant  coordinate 
algorithms  based  on  the  Lagrangian 
augmented  form  of  the  equations  of  motion 
require  the  use  of  Newton-Raphson  method 
in  order  to  ensure  that  the  constraint 
equations  are  satisfied  at  the  position  level. 
Recursive  and  joint  variables  methods  can 
also  be  used  instead  of  redundant  coordinate 
formulations  in  order  to  avoid  Newton- 
Raphson  algorithm.  Another  approach  that 
can  be  used  to  enforce  the  constraint 
equations  at  the  position  level  is  the  penalty 
method.  This  model  does  not  lead  to 


reduction  in  the  number  of  the  system 
degrees  of  freedom. 

Two  other  approaches  that  capture  the  joint 
compliance  will  also  be  considered  in  this 
study.  The  first  is  the  compliant  discrete 
element  method  that  employs  MBS  bushing 
elements  to  define  the  connectivity  between 
the  track  links.  This  approach  as  in  the  case 
of  the  penalty  method  requires  assuming 
stiffness  and  damping  coefficients  at  the 
connection,  and  therefore,  it  allows  for  the 
flexibility  between  the  track  links.  In  the 
second,  the  compliant  continuum-based  joint 
formulation  that  employs  ANCF  finite 
elements  is  used.  This  approach,  which 
captures  new  joint  deformation  modes,  leads 
to  linear  connectivity  conditions  which  can 
be  applied  at  a  preprocessing  stage  allowing 
for  an  efficient  elimination  of  the  dependent 
variables,  this  leads  to  a  constant  inertia 
matrix  and  zero  Coriolis  and  centrifugal 
forces  [18].  This  approach  leads  to  new 
types  of  FE  chain  meshes  that  have  desirable 
characteristics. 


ALGEBRAIC  CONSTRAINT 
EQUATIONS 

In  the  methods  of  constrained  dynamics , 
there  are  two  approaches  that  are  often  used 
to  model  ideal  mechanical  joints  that  do  not 
account  for  the  effect  of  elasticity  and 
damping.  These  two  methods  are  the 
augmented  formulation  that  employs  the 
technique  of  Lagrange  multipliers  or  the 
recursive  formulation  which  allows  for 
systematic  elimination  of  the  dependent 
variables  using  the  algebraic  equations. 
These  two  formulations  are  briefly  discussed 
in  this  section. 
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Augmented  Formulation 

In  the  augmented  formulation,  the  constraint 
forces  explicitly  appear  in  the  dynamic 
equations  which  are  expressed  in  terms  of 

redundant  coordinates.  The  constraint 

relationships  are  used  with  the  differential 
equations  of  motion  to  solve  for  the 
unknown  accelerations  and  constraint 

forces.  While  this  approach  leads  to  a  sparse 
matrix  structure,  it  has  the  drawback  of 
increasing  the  problem  dimensionality  and  it 
requires  more  sophisticated  numerical 

algorithms  to  solve  the  resulting  system  of 
differential  and  algebraic  equations  (DAE). 
Using  the  generalized  coordinates,  the 
equations  of  motion  of  a  body  i  can  be 
written  as  [11,  15] 


M'q' =Q:+Q:,+Q;  (1) 

where  M'is  the  mass  matrix  of  the  body, 
q'=[^R'r  ,rJ  is  the  vector  of  the 

accelerations  of  the  body  with  R'  defining 
the  body  translation  and  '  defining  the 
body  orientation,  Q[  is  the  vector  of 

external  forces,  Q  is  the  vector  of  the 
constraint  forces  which  can  be  written  in 
terms  of  Lagrange  multipliers  as 

Q'  =  -CJ  ,  C  ,  is  the  constraint  Jacobian 

c  q  q' 

matrix  associated  with  the  coordinates  of 
body  i ,  and  is  the  vector  of  the  inertia 

forces  that  absorb  terms  that  are  quadratic  in 
the  velocities.  The  constraint  equations  at 
the  acceleration  level  can  be  written  as 
CUq'  =  Q) ,  where  Q'd  is  a  vector  that 

absorbs  first  derivatives  of  the  coordinates. 
Using  equation  (1)  with  the  constraint 
equations  at  the  acceleration  level,  one 
obtains 


M  Cqr 

q 

Qe+Qv 

1 

o 

O 

1 _ 

i 

O' 

_ i 

The  matrices  and  vectors  that  appear  in  this 
equation  are  the  system  matrices  and  vectors 
that  are  obtained  by  assembling  the  body 
matrices  and  vectors.  The  preceding  matrix 
equation,  which  ensures  that  the  constraint 
equations  are  satisfied  at  the  acceleration 
level,  can  be  solved  for  the  accelerations  and 
Lagrange  multipliers.  In  order  to  ensure  that 
the  algebraic  kinematic  constraint  equations 
are  satisfied  at  the  position  and  velocity 
levels,  the  independent  accelerations  q,  are 

identified  and  integrated  forward  in  time  in 
order  to  determine  the  independent 

velocities  q,  and  independent  coordinates 
q; .  Knowing  the  independent  coordinates 
from  the  numerical  integration,  the 

dependent  coordinates  q,  can  be 

determined  from  the  nonlinear  constraint 
equations  using  an  iterative  Newton- 
Raphson  algorithm  that  requires  the  solution 
of  the  system  CqAq,;  =  -C  ,  where  Aqd  is 
the  vector  of  Newton  differences,  and  Cq 

is  the  constraint  Jacobian  matrix  associated 
with  the  dependent  coordinates.  Knowing 
the  system  coordinates  and  the  independent 
velocities,  the  dependent  velocities  qrf  can 

be  determined  by  solving  a  linear  system  of 
algebraic  equations  that  represents  the 
constraint  equations  at  the  velocity  level. 
This  linear  system  of  equations  in  the 

velocities  can  be  written  as 
cq,<5b  =-Cq,q,-Cf;  where  Cq  is  the 

constraint  Jacobian  matrix  associated  with 
the  independent  coordinates,  and 
C,  =  cC/ dt  is  the  partial  derivative  of  the 
constraint  functions  with  respect  to  time. 
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Lagrange  multipliers,  on  the  other  hand,  can 
be  used  to  determine  the  constraint  forces. 
For  a  given  joint/:,  the  generalized 
constraint  forces  acting  on  body  i , 
connected  by  this  joint,  can  be  obtained 
from  the  equation 

(q;),=-(c.)*  .=[Ff  vT  <3) 

where  F;'  and  Tlk  are  the  generalized  joint 

forces  associated,  respectively,  with  the 
translation  and  orientation  coordinates  of 
body? .  Using  the  results  of  equation  (3),  the 
reaction  forces  at  the  joint  definition  point 
can  be  determined  using  the  concept  of  the 
equipollent  systems  of  forces. 


Recursive  Formulation 


A'/ 


Figure  2:  Revolute  joint 

Another  alternate  approach  for  formulating 
the  equations  of  motion  of  constrained 
mechanical  systems  is  the  recursive  method, 
wherein  the  equations  of  motion  are 
formulated  in  terms  of  the  joint  degrees  of 
freedom.  This  formulation  leads  to  a 


minimum  set  of  differential  equations  from 
which  the  workless  constraint  forces  are 
automatically  eliminated  [11,  17].  The 

numerical  procedure  used  in  solving  these 
differential  equations  is  much  simpler  than 
the  procedure  used  in  the  solution  of  the 
mixed  system  of  differential  and  algebraic 
equations  resulting  from  the  use  of  the 
augmented  formulation.  In  the  recursive 
formulation,  the  equations  of  motion  are 
formulated  in  terms  of  joint  degrees  of 
freedom.  In  this  formulation,  the  multibody 
system  is  assumed  to  consist  of  subsystems, 
as  in  the  case  of  the  track  chains  shown  in 
figure  2  where  body  j  corresponding  to 
body  i  —  1.  The  absolute  coordinates  and 
velocities  of  an  arbitrary  body  i  in  a 
subsystem  are  expressed  in  terms  of  the 
independent  joint  variables  as  well  as  the 
absolute  coordinates  and  velocities  of  body 
i— 1.  If  body  i  is  connected  to  body  i- 1 
through  a  revolute  joint,  which  is  the  case  in 
this  subsystem,  the  relative  rotation  is  the 
only  degree  of  freedom  represented  between 
the  bodies.  The  connectivity  between  bodies 
i  and  body  i  —  lean  then  be  described  using 
the  kinematic  relationships 

R'  +  A'u'n  -  R'  1  -  A'“'u'p 1  =  o] 

}  (4) 

where  R'  is  the  global  position  vector  of  the 
origin  of  body  i;  A'  is  the  transformation 
matrix  that  defines  the  body  orientation  and 
can  be  expressed  in  terms  of  Euler 
parameters;  uj,  and  uj. 1  are  the  local 

position  vectors  of  point  P  defined  in  the 
coordinate  systems  of  body  i  and  i-1, 
respectively,  ‘  and  1-1  are,  respectively, 
the  absolute  angular  velocity  vectors  of 


bodies  i  and  /  —  1,  and  ,  '~1  is  the  angular 
velocity  vector  of  body  i  with  respect  to 
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body  i- 1  which  can  be  defined  as 
'■'-1  =(|)'v'_1  with  v'_1  =  v'_1  where 

t 

v'1  is  a  unit  vector  along  the  axis  of 
rotation  defined  in  the  coordinate  system  of 
body  i— 1,  and  <()'  is  the  angle  of  relative 
rotation.  By  differentiating  the  first  equation 
in  equation  (4)  twice  and  the  second  once 
with  respect  to  time,  one  obtains 

R  =-  !  X(  1  XUp)  —  (  '  X Up) 

+  Rm+  mx(  xu';')  +  (  Hxu'-')  (5) 

*'=  M  +  vMf+(  MxvM)f 

In  this  equation,  k  is  the  absolute  angular 
acceleration  vector  of  body  k  .  Using  the 
kinematic  equations  obtained  in  this  section, 
one  can  systematically  eliminate  the 
dependent  variables  in  order  to  obtain  a 
number  of  differential  equations  of  motion 
equal  to  the  number  of  the  system  degrees  of 
freedom.  Using  this  approach,  one  obtains  a 
dense  inertia  matrix  in  a  system  of  dynamic 
equations  that  does  not  have  constraint 
forces.  A  second,  alternative  approach  is  to 
use  the  kinematic  equations  developed  in 
this  section  to  determine  all  the  system 
absolute  coordinates  and  velocities.  One  can 
then  construct  equation  (2),  which  can  be 
solved  for  the  accelerations  and  Lagrange 
multipliers.  Using  the  absolute  acceleration 
relationships  of  equation  (5),  one  can 
determine  the  relative  joint  accelerations. 
The  joint  accelerations  can  be  integrated 
forward  in  time  in  order  to  determine  the 
joint  coordinates  and  velocities. 


generalized  and  the  Cartesian  moments  [11, 
17].  This  is  important  in  interpreting  the 
reaction  forces  of  the  ideal  joints  and  also 
important  in  the  implementation  of  the 
penalty  method  and  bushing  elements.  Let 
F'  be  a  force  vector  that  acts  at  a  point  P‘ 
on  a  rigid  body  i.  If  this  force  vector  is 
assumed  to  be  defined  in  the  global 
coordinate  system,  then  the  virtual  work  of 
this  force  vector  can  be  written  as 

§Wle  =  F,78rj,,  where  8rj,  can  be  found 
using  the  virtual  change  in  the  position 
vector  of  an  arbitrary  point  on  rigid  body  i 
as 

—  r  8  r  ' 

8r;=[l  -A'OpG']  .  (6) 

In  this  equation,  A' is  the  transformation 
matrix  that  defines  the  body  orientation,  uj, 
is  the  skew  symmetric  matrix  associated 
with  the  vector  Up  that  defines  the  local 

coordinates  of  the  point  Pl ,  and  G‘  is  the 
matrix  that  relates  the  angular  velocity 

vector  1  defined  in  the  body  coordinate 
system  to  the  time  derivatives  of  the 

orientation  coordinates,  that  is  '  =G'  ' . 

_  •  ■  •  -T 

Note  that  since  Up  =  A'Up  A'  ,  equation  (6) 

can  be  written  as  8rj,  =SR'  -UpG'S  ' 

Using  this  equation  in  the  virtual  work 
expression,  one  obtains 

Sir;  =F'78R' -Fi'u'/)G'8  ',  which  can  be 
written  as 


GENERALIZED  FORCES 

In  defining  the  joint  forces  between  the  track 
links,  it  is  important  to  understand  the 
relationship  and  differences  between  the 

Page 


8w;=f;sr' +f'7s  ‘  (7) 

where  Fj.  =  F' ,  and  Fg=-G'  u'pF'.  These 
equations  imply  that  a  force  that  acts  at  an 
arbitrary  point  on  the  rigid  body  i  is 
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equipollent  to  another  system  defined  at  the 
reference  point  that  consists  of  the  same 
force  and  a  set  of  generalized  forces,  defined 

by  Fg=-G'u'pF'  associated  with  the 
orientation  coordinates  of  the  body  [11,  17]. 

Since  u'p  is  a  skew-symmetric  matrix,  it 

■T 

follows  that  u'p  =  -u'p  .  Using  this  identity, 
one  can  show  that  the  generalized  moment 
can  be  written  as  Fg  =G'  (u'pxF')  ,  where 
IVF  =UpxF'  is  the  Cartesian  moment 
resulting  from  the  application  of  the  force 
F',  and  G'  is  the  matrix  that  relates  the 
angular  velocity  vector  '  defined  in  the 
global  coordinate  system  to  the  time 
derivatives  of  the  orientation  coordinates, 

that  is  '  =  G'  ' .  It  follows  that  the 
relationship  between  the  generalized  and 

Cartesian  moment  is  Fg  =G'  If  the 
components  of  the  moment  vector  are 
defined  in  the  body  coordinate  system,  one 

has  Fg  =G''m',  where  M'  =A''M^ 

Relationships  developed  in  this  section  will 
be  used  in  the  formulation  of  the  joint  forces 
in  the  case  of  the  penalty  method.  These 
relationships  will  also  be  used  in  the 
computer  implementation  of  the  bushing 
element  in  general  MBS  algorithms. 


JOINT  FORMULATIONS 

In  this  paper,  four  different  joint  models  are 
considered;  two  models  are  based  on 
algebraic  equations  that  require  the  use  of 
the  methods  of  constrained  dynamics  or  the 
penalty  method.  In  the  case  of  constrained 
dynamics,  an  alternative  to  the  use  of 
Lagrange  multipliers  is  the  use  of  the 


recursive  methods,  as  previously  discussed. 
When  Lagrange  multiplier  technique  or  the 
recursive  methods  are  used  the  constraint 
equations  must  be  satisfied  at  the  position, 
velocity,  and  acceleration  levels.  The 
penalty  method,  on  the  other  hand,  cannot 
satisfy  the  algebraic  constraint  equations  at 
the  acceleration  level. 


Constraint  Equations 

The  revolute  (pin)  joint  is  used  in  this 
section  as  an  example  to  demonstrate  the 
formulation  of  the  algebraic  constraint 
equations.  This  joint  has  been  used  in  the 
literature  in  the  modeling  of  the  track  chains. 
As  shown  in  figure  2,  the  track  chain  can  be 
assumed  to  consist  of  links  connected  to 
each  other  by  a  pin  joint  that  allows  for  the 
relative  rotation  between  them.  In  general 
MBS  algorithms,  the  nonlinear  algebraic 
equations  that  define  the  pin  joint  are 
expressed  in  terms  of  the  absolute 
coordinates  of  the  two  bodies  i  and  j 
connected  by  the  joint.  The  five  algebraic 
constraint  equations  that  eliminate  five 
degrees  of  freedom  can  be  written  in  terms 
of  the  absolute  Cartesian  coordinates  of  the 
two  bodies  as  [17] 


C(q\qJ)- 


|v,  v-1 


v'2v' 


V  r 

M  A/ 


V  rlJ 

v2  lP 


=0 

(8) 


where  v'  and  V  are  two  vectors  defined 
along  the  joint  axis  on  bodies  i  and  j, 

respectively;  v',v],v'2  form  an  orthogonal 
triad  defined  on  body  i ;  Ag  is  a  constant 
[17];  and 


v'p  =  y'p-yp  =R'  +A'Up-R'  -A'Up  (9) 
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In  this  equation,  A'  and  Aj  are  the 
transformation  matrices  that  define  the 
orientation  of  bodies  i  and  j ,  respectively, 
and  u‘P  and  uJp  are  the  local  position  vectors 
of  points  P'  and  Pj  with  respect  to  bodies 
i  and  j ,  respectively.  Points  P'  and  Pj  are 
defined  on  the  axis  of  the  pin  joint  on  bodies 
i  and  j ,  respectively.  One  can  show  that 
the  Jacobian  matrix  of  the  pin  joint 
constraints  is  defined  as 


C 


<j 


v,7h; 

v'7H' 

v'7H( 

-fHj  +  v'fHj, 

-v'fHy 

fH'+v'lH' 

-v'THy 

2ifH'P 

(10) 


where  C ,  and  C  ,  are  the  constraint 

q  q 

Jacobian  matrices  associated  with  the 
coordinates  of  bodies  i  and  j ,  respectively; 
and  other  vectors  and  matrices  that  appear  in 
the  preceding  equation  are 


H'd  =  [I 


Air;  G'],  Hj,  =  [i 


AJuJP  Gy], 


H;=5y'  =  d  (Ay\),  H,  =  5v2=A (AV2), 
1  dq  dq1  1  2  cq‘  ctf 


dq1  dq1 

Another  alternate  approach  for  formulating 
the  revolute  joint  constraints  is  to  consider  it 
as  a  special  case  of  the  spherical  joint  in 
which  the  relative  rotation  between  the  two 
bodies  is  allowed  only  along  the  joint  axis. 
If  point  P  is  the  joint  definition  point,  and 
v'  and  v'  are  two  vectors  defined  along  the 
joint  axis  on  bodies  i  and  j  ,  respectively, 


the  constraint  equations  of  the  revolute  joint 
can  be  written  as 

C(q',qy)  =  [rp  v'[vJ  v'jV7]  =0.  The 

last  two  equations  in  this  equation  guarantee 

that  the  two  vectors  v'  and  \J  remain 
parallel,  thereby  eliminating  the  freedom  of 
the  relative  rotation  between  the  two  bodies 
in  two  perpendicular  directions. 


Penalty  Method 

The  constraint  equations  that  describe  the 
connectivity  between  the  track  links  can  be 
enforced  using  the  penalty  approach.  In  this 
case,  these  algebraic  equations  are  not 
satisfied  at  the  acceleration  level.  The 
penalty  method  does  not  lead  to  elimination 
of  degrees  of  freedom,  and  therefore,  it  is 
conceptually  different  from  the  case  of 
Lagrange  multiplier  technique  or  the 
recursive  approach.  In  order  to  demonstrate 
the  penalty  approach,  the  violations  in  the 
constraint  equations  of  a  revolute  joint  k 
can  be  written  as 


d«r=[rP  VTV2  V* ( V^  ]  (11) 

Using  this  violation  dk  ,  a  restoring  force 
vector  can  then  be  defined  as 
fk=hkdk+ckdk,  where  kk,  and  ck  are 
assumed  penalty  stiffness  and  damping 
coefficients,  respectively,  and  dk  is  the  time 
derivative  of  the  violation  vector  d^ .  The 
virtual  work  of  this  restoring  force  fk  can 
then  be  written  as  8Wk  =-fk‘bdk,  which 
can  be  written  as 
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8 W?  =  -Fj8r«  -  Fxh  (vf  v7 )  -  F2 8  (v2  v7 ) 

(12) 


proper  selection  of  the  penalty  coefficients, 
these  forces  ensure  that  the  constraint 
equations  are  satisfied  at  the  position  level. 


where 

8rj  =  8R  -u'BG'S  '-8R7+u7G78  7, 

Fb  =  kkr‘i  +  ckr‘i ,  8(v|V)  =  v7T8v|  +v|r8v;' , 
8(\'2t\j)  =  v,ySv'2  +  v27Sv7 , 
d(\'Tx\J) 


+  C u 


f.=*.Kv') 


dt 

j(v'2vJ) 


,  and 


dt 


Equation  (12)  can  be  used  to  define  a  set  of 
generalized  forces  acting  on  bodies  i  and  j 
that  maintain  the  connectivity  between  the 
two  bodies  as 

8  W?  =  Q*SR  +  Q;' 89 '  +  Q(8R '  +  80 j 

(13) 


Which  can  be  rewritten  in  a  compact  form  as 
8W=Q-8q'+Qir8q^  where  q  and  q' 

are  the  generalized  coordinates  of  bodies  1 
and  J  ,  and 


Qr 

-F* 

_Qe_ 

g,tu7fb+g'tm;+g,tm2 

oi 

F 

-GjTuiTFB  -g,7m;  -G'tM, 

(14) 


where  Mj  =  -F1v|v' ,  M/  =  Fx\ |  vJ , 
M'2  =  -FjV'2  v' ,  and  M2  =  Ejv2  \J . 

Equation  (14)  defines  the  generalized  forces 
associated  with  the  absolute  Cartesian 
coordinates  due  to  the  revolute  joint 
connection  between  bodies  i  and  j .  With  a 


COMPLIANT  DISCRETE  ELEMENT 
JOINT  FORMULATION 


Figure  3:  Bushing  element 

The  compliant  discrete  element  joint 
formulation  allows  for  introducing  joint 
deformations.  In  this  approach,  no  algebraic 
equations  are  enforced.  In  most  MBS 
computer  codes,  the  compliant  discrete 
element  joint  method  can  be  applied  using 
the  standard  MBS  bushing  element  that 
allows  for  introducing  three  force  and  three 
moment  components  that  can  be  linear  or 
nonlinear  functions  of  the  body  coordinates. 
As  shown  in  figure  3,  the  position  vectors 

u/,  and  Up  of  two  points,  Px  and  P2J  on 

body  j,  can  be  used  to  define  one  axis  of  the 
coordinate  system  of  the  bushing  element  as 


ny 


sT 

1 

Cl 

u7  -u7 

V  Pi  Pi  // 

P\  Pi 

where  n 1  is 


one  of  the  bushing  axes  defined  in  the  body 
j  coordinate  system.  This  axis  can  then  be 


Page  12  of  22 


CHAIN  DYNAMIC  FORMULATIONS  FOR  MULTIBODY  SYSTEM  TRACKED  VEHICLES 


UNCLASSIFIED 


Proceedings  of  the  2012  Ground  Vehicle  Systems  Engineering  and  Technology  Symposium  (GVSETS) 


defined  in  the  global  coordinate  system  as 
n  =  A7n7,  where  A7  is  the  transformation 
matrix  that  defines  the  orientation  of  the 
coordinate  system  of  body  j  in  the  global 
system.  Using  this  axis,  one  can  define  the 
directional  properties  of  the  bushing  element 
with  the  other  two  axes  of  the  bushing 
coordinate  system  being  defined  using  the 

transformation  matrix  Abj  =[t/  t/  n 

where  t{  and  t7  are  the  two  unit  vectors 
that  complete  the  three  orthogonal  axes  of 
the  bushing  element  coordinate  system. 
Assuming  that  body  j  is  a  rigid  body,  the 
bushing  coordinate  system  can  be  defined 
with  respect  to  the  global  coordinate  system 
as  Abj  =  A 7  A hj .  Choosing  points  P'  and 
Pj7  to  initially  coincide;  one  can  define  the 
bushing  deformation  and  rate  of  deformation 
vectors  in  the  bushing  coordinate  system  as 

~bii=A bjTr\  and  ' Mj -- Ahy  r,J , 

respectively,  where  r'7  =  r'p  -  r7  is  the 

position  vector  of  point  P1 7  with  respect  to 
point  P‘ . 

The  rotational  deformation  of  the  bushing 
element  can  be  obtained  using  the 
transformation  matrix  that  defines  the 
orientation  of  the  bushing  coordinate  system 
on  body  i  with  respect  to  the  bushing 
coordinate  system  on  body  j  .  This  matrix  is 

defined  as  A hij  =  Ahj  A1”  where  Abj  is  the 
orientation  matrix  of  the  bushing  coordinate 
system  on  body  j,  while  Abl  is  the 
orientation  matrix  of  the  bushing  coordinate 
system  with  respect  to  the  coordinate  system 
of  body  i  that  is  defined  as  Abl  =  A‘Ab' . 
Assuming  that  the  relative  rotations  between 
bodies  i  and  j  are  small,  the  relative 

rotation  matrix,  Ahl]  can  be  used  to  extract 


three  relative  rotations  defined  in  the 
bushing  coordinate  system, 

~bij  =[  f  byij  fj  The  relative 

angular  velocity  between  the  two  bodies 
defined  in  the  bushing  coordinate  system 

can  also  be  written  as  ~bii  =  Ab'  (  '  -  7  ), 

where  '  and  7  are  the  absolute  angular 
velocity  vectors  of  bodies  i  and  j , 
respectively,  defined  in  the  global 
coordinate  system.  The  bushing  stiffness  and 
damping  coefficients  are  often  determined 
using  experimental  testing,  and  these 
coefficients  are  defined  generally  in  the 
bushing  coordinate  system.  Let  K  r  and  C  r 
be  the  translational  stiffness  and  damping 
matrices,  respectively,  defined  with  respect 
to  the  bushing  coordinate  system;  and 
assume  that  the  rotational  stiffness  and 
damping  matrices  are  Ke  and  C0 , 

respectively.  In  terms  of  translational  and 
rotational  stiffness  and  damping  matrices, 
the  force  vector  defined  in  the  bushing 
coordinate  system  can  be  written  as 


X 1 

X 

o " 

bij 

+ 

X 

o " 

bij 

= 

_Me  _ 

0 

* 

<x> 

1 _ 

bij 

0 

— i 

<x> 

V 

—  bij 

(15) 


This  force  vector  can  then  be  defined  in  the 
global  coordinate  system,  and  the  results  can 
be  used  to  define  the  generalized  bushing 
forces  and  moments  acting  on  the  two 
bodies  as  previously  described  in  this  paper. 
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COMPLIANT  CONTINUUM-BASED 
JOINT  FORMULATION 


Figure  4:  ANCF  beam  element  coordinates 

The  compliant  continuum-based  joint 
formulation  allows  capturing  joint  strain 
modes  that  cannot  be  captured  using  the 
compliant  discrete  element  joint  method. 
When  ANCF  finite  elements  are  used,  one 
can  develop  new  FE  meshes  that  have  linear 
connectivity  and  constant  inertia  [18].  This 
allows  for  systematically  eliminating 
dependent  variables  at  a  preprocessing  stage, 
and  as  a  result,  there  is  no  need  for  the  use 
of  joint  formulations  in  the  main  processor. 
This  approach  can  be  used  to  develop  new 
spatial  chain  models  where  the  modes  of 
deformations  at  the  definition  points  of  the 
joints  that  allow  for  rigid  body  rotations 
between  ANCF  finite  elements  can  be 
captured.  The  displacement  field  of  an 
ANCF  finite  element,  as  the  one  shown  in 
figure  4,  can  be  written  as 
r(x,y,z,t)=s(x,y,z)e(t)  where  x,y,  and 
Z  are  the  element  spatial  coordinates;  t  is 
time;  S  is  the  element  shape  function 
matrix,  and  e  is  the  vector  of  element  nodal 
coordinates.  Using  this  displacement  field, 


the  equations  of  a  pin  joint  between 
elements  i  and  j  can  be  written  using  the 

six  scalar  equations  r' =  r7,  =  rj  ,  where 

a  is  the  coordinate  line  that  defines  the 
joint  axis;  a  can  be  x,  y,  or  z  or  any  other 
coordinate  line.  The  six  scalar  equations 
eliminate  six  degrees  of  freedom;  three 
translations,  two  rotations,  and  one 
deformation  mode.  Therefore,  this  joint  has 
five  modes  of  deformation  that  include 
stretch  and  shear  modes.  This  ANCF 

revolute  joint  model  ensures  C1  continuity 
with  respect  to  the  coordinate  line  a  and 
C°  continuity  with  respect  to  the  other  two 
parameters.  It  follows  that  the  Lagrangian 

strain  component  £aa  =  (rjra -l)/2  is 

continuous  at  the  joint  definition  point, 
while  the  other  five  strain  components  can 
be  discontinuous.  The  resulting  joint 
constraint  equations  are  linear,  and 
therefore,  can  be  applied  at  a  preprocessing 
stage  to  systematically  eliminate  the 
dependent  variables.  Using  these  equations, 
one  can  develop  a  new  kinematically  linear 
FE  mesh  for  flexible-link  chains  in  which 
the  links  can  have  arbitrarily  large  relative 
rotations. 

The  comparative  numerical  study  presented 
in  the  following  section  will  be  focused  on 
three  models;  the  ideal  constrained  joint 
model,  the  ideal  penalty  joint  model,  and  the 
compliant  discrete  element  model  (bushing). 
The  results  of  the  ANCF  joint  model  will  be 
considered  in  future  investigations  in  order 
to  allow  for  presenting  a  more  meaningful 
study  that  compares  different  ANCF  joint 
formulations. 
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SIMULATION  RESULTS 


Figure  5:  Tracked  vehicle  different  contact 
models 


Figure  6:  Suspension  system  layout  of 
Ml  13  tracked  vehicle  [6] 

In  this  section,  numerical  results,  obtained 
using  Ml  13  tracked  vehicle  model  shown  in 
figure  1,  are  used  to  compare  the  different 
joint  formulations  presented  in  this 
investigation.  The  Ml  13  tracked  vehicle  is 
armored  personnel  carrier  that  consists  of  a 
chassis,  idler,  sprocket,  5  road-wheels,  and 
64  track  links  on  each  track  side  (right  and 
left).  Figure  5  shows  the  engagement  of  the 
track  links  with  some  of  the  vehicle 


components,  while  more  details  about  the 
roller  arrangement  and  the  configuration  of 
the  suspension  system  of  the  Ml  13  tracked 
vehicle  model  are  shown  in  figure  6.  Table  1 
shows  the  inertia  prosperities  for  all  the 
tracked  vehicle  model  components  used  in 
this  simulation.  More  specifications  of  this 
vehicle  can  be  obtained  from  open  sources 
[5],  The  vehicle  has  a  suspension  system 
that  consists  of  road  arms  placed  between 
the  road  wheels  and  chassis  as  well  as  shock 
absorbers  connected  to  each  road  arm.  Table 
2  shows  the  stiffness  and  the  damping 
coefficients  of  the  contact  models  as  well  as 
their  friction  coefficients. 

In  this  study,  two  different  simulation 
scenarios,  one  with  the  suspension  system 
and  the  other  without  the  suspension  system, 
will  be  considered  to  study  the  effect  of 
using  the  suspension  system  on  the  results. 
The  road  arms  and  the  sprockets  are 
connected  to  the  chassis  by  revolute  joints, 
and  the  road  arms  are  connected  to  the  road- 
wheels  by  revolute  joints.  The  track  links  are 
connected  to  each  other  using  re  volute 
joints,  which  can  be  modeled  using  the 
penalty  method,  bushing  element,  or  ANCF 
finite  elements  as  previously  mentioned. 
Tensioners  are  added  to  the  system  by 
connecting  each  idler  to  a  tensioner  with  a 
re  volute  joint  and  connecting  the  tensioner 
to  the  chassis  with  a  prismatic  joint  to  ensure 
only  translation  between  them. 
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Figure  7:  Sprocket  angular  velocity 

The  angular  velocities  of  the  sprockets  of 
the  vehicle  model  considered  in  this 
numerical  investigation  are  assumed  to 
increase  linearly  after  1  sec  until  it  reaches  a 
constant  value  of  25  rad/sec  after  8  seconds 
as  shown  in  figure  7.  Figure  8  shows  the 
chassis  vertical  displacement  using  the  joint 
model  with  and  without  the  suspension 
system.  The  results  presented  in  this  figure 
show  that  the  model  with  a  suspension 
system  has  less  vibration  compared  to  the 
model  without  suspension. 

The  penalty  method  model  and  bushing 
element  model  shown  in  figures  9-14  each 
have  a  stiffness  coefficient  of  109  N/m.  In 
real  life  applications,  the  stiffness  coefficient 
of  bushing  elements  used  in  track  chain 
connections  are  much  lower;  they  are 
increased  in  this  case  for  better  comparisons 
to  the  re  volute  joint  and  penalty  models.  The 
rotational  stiffness  and  damping  coefficients 
about  the  axis  of  rotation  for  the  track  links 
of  the  bushing  element  model  are  set  to  zero 
in  this  case  to  allow  for  rotation  about  the 
lateral  axis.  Figures  9  and  10  show, 
respectively,  the  chassis  forward  position 


and  velocity  results  obtained  using  different 
joint  models.  While  the  results  presented  in 
these  figures  show  good  agreement,  the 
computational  time  varies  when  these 
different  models  are  used.  The  constrained 
joint  model  takes  less  computational  time 
compared  to  the  penalty  and  the  bushing 
element  models;  the  penalty  model  CPU 
time  is  twice  that  of  the  constraint  model 
CPU  time,  while  the  bushing  model  CPU 
time  is  four  times  that  of  the  constraint 
model  CPU  time.  This  increase  in  the 
bushing  model  CPU  time  is  attributed  to  the 
high  stiffness  coefficients  used  in  this 
model. 


(—  without  suspension,  with 
suspension) 
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Time  (s) 

Figure  9:  Chassis  forward  position 


—  Constrained  joint  model,  — ▼ — 
method  model,  -fi-  Bushing  element 
model) 


Figure  10:  Chassis  forward  velocity 

(  A  Constrained  joint  model,  — ▼ — 
Penalty  method  model,  — 0-- Bushing  element 
model) 


Figure  1 1  shows  the  vertical  displacement  of 
a  track  link  in  the  global  coordinate  system 
using  different  joint  models.  While  the 
results  show  good  agreement  the  joint  model 
has  less  vibration  compared  to  the  penalty 
and  bushing  element  models.  Figure  12 
shows  the  motion  trajectory  of  a  track  link  in 
the  chassis  coordinate  system  using  the 
constrained  model,  the  penalty  model,  and 
the  bushing  element  model,  respectively. 
While  the  results  presented  in  these  figures 
show  good  agreement,  it  is  clear  that  the 
constrained  dynamics  model  solution 
exhibits  less  oscillations. 


Figure  11:  Vertical  displacement  of  a  track 
link 


(  A  Constrained  joint  model,  — ▼ — 
Penalty  method  model,  — Bushing  element 
model) 
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Local  Y-Position  (Longitudinal,  m) 


Figure  12:  Trajectory  motion  of  a  track  link 
in  the  chassis  coordinate  system 

(  Constrained  joint  model,  Penalty 
method  model,  Bushing  element  model) 
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Figure  13:  Joint  longitudinal  forces 

(  Constrained  joint  model,  Penalty 
method  model,  k  =  lxlO9  N/m) 


Figures  13  -  16  show  the  joint  forces  in  the 
longitudinal  and  the  vertical  directions, 
respectively,  obtained  using  the  constrained 
and  penalty  joint  models.  The  penalty  results 
are  obtained  using  a  stiffness  coefficient  of 

lxlO9  N/m  and  damping  coefficient  of 
5xl04  N.s/m.  The  results  show  good 
agreement  between  the  constrained  and 
penalty  joint  models.  Figures  15  and  16 
show  the  same  results  in  the  case  of  a 
stiffness  coefficient  of  lxlO7  N/m  and 
damping  coefficient  of  5xl04  N.s/m. 


Figure  14:  Joint  vertical  forces 

(  Constrained  joint  model,  —  Penalty 
method  model,  k  =  lxlO9  N/m) 
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Figure  15:  Joint  Longitudinal  forces 


(  Constrained  joint  model,  —  Penalty 
method  model,  k  =  lxlO7  N/m) 


Figure  16:  Joint  vertical  forces 

(  Constrained  joint  model,  —  Penalty 
method  model,  k  =  lxlO7  N/m) 


The  results  of  figures  15  and  16,  which 
show  significant  differences  between  the 
two  models,  demonstrate  the  drawback  of 
the  penalty  method  when  the  penalty 
stiffness  coefficient  is  reduced.  Similar 
results  can  be  expected  in  the  case  of  the 
bushing  element  models  where  the  forces 
obtained  also  depend  on  the  stiffness  and 
damping  coefficient  of  the  joint.  Figure  17 
shows  the  joint  deformation  predicted  using 
the  penalty  model  using  different  stiffness 
coefficients.  The  penalty  model  with  a 
stiffness  coefficient  of  lxlO9  N/m  shows 
much  less  deformation,  less  than  0.06  mm, 
while  the  penalty  model  with  a  stiffness  of 
lxlO7  N/m  has  much  more  deformation, 
over  2.5  mm,  between  the  track  links. 


Figure  17:  Joint  deformation  using  penalty 
model 


(  k  =  lxlO9  N/m,  —  k  =  lxlO7  N/m) 
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SUMMARY  AND  CONCLUSIONS 

In  this  investigation,  different  MBS  joint 
formulations  are  presented  and  compared 
using  detailed  tracked  vehicle  models.  Three 
main  joint  formulations  are  discussed;  they 
are  the  ideal  joint  formulation,  the  compliant 
discrete  element  joint  formulation,  and  the 
compliant  continuum-based  joint 
formulation.  The  ideal  joint  formulation  is 
developed  to  eliminate  the  relative 
displacement  between  the  two  bodies 
connected  by  the  joint.  This  can  be  achieved 
by  enforcing  a  set  of  joint  algebraic 
equations  using  a  constrained  dynamics 
approach  or  by  using  the  penalty  method. 
The  constrained  dynamics  approach 
eliminates  degrees  of  freedom  and  ensures 
that  the  constraint  equations  are  satisfied  at 
the  position,  velocity,  and  acceleration 
levels.  The  penalty  method,  on  the  other 
hand,  does  not  reduce  the  number  of  degrees 
of  freedom  and  ensures  that  the  constraint 
equations  are  satisfied  at  the  position  level 
only  provided  that  a  high  stiffness 
coefficient  is  used.  The  compliant  discrete 
element  formulation,  which  allows  for  joint 
deformations,  can  be  systematically  applied 
using  a  standard  MBS  bushing  element  that 
allows  for  six  degrees  of  freedom  of  relative 
motion.  The  compliant  continuum-based 
approach  can  be  used  to  develop  new  joints 
that  capture  deformation  modes  that  are  not 
captured  by  the  compliant  discrete  element 
joint  formulation.  ANCF  finite  elements  can 
be  used  to  systematically  develop  new  joints 
with  distributed  elasticity  and  linear 
connectivity  conditions.  As  discussed  in  this 
paper,  it  is  important  to  choose  the  proper 
stiffness  and  damping  coefficients  when  the 
penalty  method  and  the  bushing  elements 
are  used. 


Numerical  results  were  presented  in  order  to 
compare  between  different  methods.  The 
ideal  joint  formulation  produces  the  desired 
joint  kinematics  and  accurate  joint  forces. 
The  same  is  true  with  the  penalty  force 
based  joint  when  large  penalty  stiffnesses 
are  used.  Penalty  force  based  joint 
construction  has  been  shown  to  be  sensitive 
to  the  selection  of  penalty  stiffness  with  the 
higher  stiffness  coefficients  leading  to  better 
overall  results.  However,  higher  penalty 
stiffness  increases  CPU  time  significantly 
due  to  higher  frequencies.  The  penalty 
method  and  bushing  element  models  each 
have  much  larger  CPU  times  than  the  ideal 
constrained  model  due  to  these  high  stiffness 
coefficients.  Bushing  force  based  joint 
formulation  when  used  with  very  large 
stiffnesses  has  been  shown  to  not 
automatically  equate  to  the  penalty  force 
based  joint  model.  The  ANCF  joint 
formulation  is  being  implemented  and  will 
be  compared  with  the  previous  models  in 
future  work. 
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Table  1:  Mass  and  inertia  values  for  tracked  vehicle  parts 


Part 

Mass  (kg) 

Ixx  ( kg- nr) 

Iyy  ( kg. nr) 

Izz  ( kg. nr) 

Ixy,  Iyz,  Ixz  (kg. nr) 

Chassis 

5489.2 

1786.9 

10450 

10721 

0 

Sprocket 

436.67 

13.868 

12.216 

12.216 

0 

Idler 

429.57 

14.698 

12.545 

12.545 

0 

Road  Wheel 

561.07 

26.063 

19.819 

19.819 

0 

Road  Arm 

28.016 

0.30241 

0.16252 

0.18809 

0 

Track  Fink 

39.948 

0.091164 

0.49787 

0.55977 

0 

Table  2:  Contact  Parameters 


Parameters 

Sprocket-Track  Contact 

Roller-Track  Contact 

Ground-Track  Contact 

k 

2.00xl06  N/m 

2.00xl06  N/m 

2.00xl06  N/m 

c 

5.00x10"  N  s/m 

5.00xl03  N  s/m 

5.00x10"  Ns/m 

0.150 

0.100 

0.300 
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